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The critical behavior of collective modes and the collapsing 
dynamics of trapped Bose-Einstein condensates with attrac- 
tive interactions are studied analytically and numerically. The 
time scales of these dynamics both below and above the crit- 
ical point of the collapse are found to obey power laws with 
a single parameter of N/N c — 1, where N is the number of 
condensate atoms and N c is the critical number. The collaps- 
ing condensate eventually undergoes rapid implosion, which 
occurs several times intermittently, and then the implosion 
turns to an explosion. The release energy of the explosion 
is found to be proportional to the square of the interaction 
strength, inversely proportional to the three-body recombina- 
tion rate, and independent of the number of condensate atoms 
and the trap frequency. 
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I. INTRODUCTION 

Bose-Einstein condensates (BECs) of trapped atomic 
vapor have been realized in several species. The static 
and dynamical properties of BEC crucially depend on 
the sign of the interatomic interaction, which is, with- 
out an external field, repulsive for 87 Rb |Q, 23 Na [pi, and 
*H [§, and attractive for 7 Li [g and 85 Rb §]. The in- 
teractions can be manipulated by applying the electric 
or magnetic field, as the scattering length is sensitive to 
the external field near the Feshbach resonance || . Using 
this technique, one can control not only the strength of 
the interaction ]7| but also its sign 

In a spatially uniform three-dimensional system, the 
attractive interaction between atoms makes BEC unsta- 
ble. In a spatially confined system, however, the zero- 
point energy serves as a kinetic obstacle against collapse, 
allowing metastablc BEC to be formed if the number N 
of BEC atoms is below a certain critical value N c . 
Various properties of BEC below and above N c have been 
predicted. Just below N c , BEC may colla pse via macro- 
scopic quantum tunneling 



15 



When N ex- 



ceeds N c , BEC becomes unstable and exhibits compli- 
cated collapsing dynamics. Kagan et al. |2l| have pre- 
dicted that the implosion of BEC eventually turns to 
explosion due to loss of atoms by three-body molecu- 
lar recombination. Because the lost atoms are replen- 
ished by a thermal cloud, the number of BEC atoms 
again exceeds N c , leading to the collapse- and-growth cy- 



cles of BEC [ p2p^j2l| . Such dynamic behavior of col- 
lapsing BEC has been explored with recent Rice exper- 
iments p3]. Ueda and Huang |2j| have shown that the 
collapse occurs locally within a small region around the 
center of the trap, and this result has been verified nu- 
merically [^(|. We have predicted that implosions 
occur intermittently and in rapid sequence due to com- 
petition between a loss of atoms and the attractive inter- 
actions. When the interaction is suddenly switched from 
repulsive to attractive by using the Feshbach resonance, 
various patterns such as a shell structure may be formed 
in the atomic density [^7) . 

In recent experiments with 85 Rb at JILA H, nearly 
pure condensates have been produced in which the num- 
ber of atoms can be fixed, and both the strength and sign 
of the interaction can be controlled by using the Feshbach 
resonance. These results suggest that N c can be tuned 
with the number of BEC atoms held fixed. We may thus 
study the critical behavior of BEC by controlling the rel- 
evant parameters. 

In the present paper, we study the dynamics of BEC 
with attractive interactions at zero temperature both be- 
low and above the critical point. We find that both 
collective-mode frequencies just below the critical point 
and the time it takes BEC just above this critical point 
to collapse obey power laws with a single parameter 
N/N c — 1. We use a time-dependent Gaussian trial wave 
function to analytically study collective-mode frequen- 
cies and the collapsing dynamics, and compare the ob- 
tained results with those obtained by numerically solving 
the Bogoliubov equations and the time-dependent Gross- 
Pitaevskii (GP) equation ^8|. When the atomic loss due 
to inelastic collisions is taken into account, the collapsing 
condensate undergoes intermittent implosions followed 
by an explosion. We find the energy released by the ex- 
plosion is proportional to g 2 /L^ and independent of the 
number of atoms and the trap frequency, where g is the 
strength of the interaction and L3 is a coefficient of the 
three-body recombination rate. 

This paper is organized as follows. Section || formu- 
lates the problem and introduces an effective-potential 
model based on a time-dependent Gaussian variational 
wave function. Section III discusses the collective oscil- 



lations of BEC below N c . Section IV analyzes the col- 
lapsing dynamics of BEC above N c . In both regimes, our 
analytical and numerically exact results are in excellent 
agreement. Section [y] discusses the dynamics of implo- 
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sion and explosion of BEC by incorporating the atomic 



loss into our theory. Finally, section VI summarizes our 
results. 



II. FORMULATION OF THE PROBLEM 

We consider a system of Bose-condensed atoms with 
mass m and s-wave scattering length a, confined in an 
isotropic parabolic potential. The Hamiltonian for the 
system is given by 



H 



dr ffi(r) 



2m 



27rh 2 a 



"ijr (r)^(r) 



(1) 



where ^(r) is the field operator of the atoms. The transi- 
tion amplitude of the system from an initial ipi to a final 
state ipf is expressed in terms of path integrals as 

(^fle^^-^lVi) = J Vi\)V^*e* s ^' r \ (2) 

where the action S[ip, ip*] is given by 

S[i>,ip*} =Nh Jdrjdt 



Here the length, time, and ifj are normalized in units of 
d = (h/mtoo) 1 / 2 , uJq 1 and (iV/dp) 1 / 2 , respectively, and 
g = AirNa/do. The wave function is then normalized 
to unity. Separating ip(r,t) into amplitude and phase as 
ip(r, t) = A(r, t)e l ^ r ^ , and substituting this into Eq. (g), 
we obtain 



S = Nh dr dt 



iAA - A 2 4> + ^AV 2 A - ^A 2 {V<t>) 2 



- l -V(A 2 ^)- r -A 2 - 9 -A± 



(4) 



The requirement that the action be stationary with re- 
spect to small variations in (f> leads to 



OA 2 



V(A 2 V0) = 0. 



(5) 



This is nothing but the equation of continuity that guar- 
antees the conservation of the number of atoms. The 
stationary condition for small variations in A is, together 
with Eq. (||), equivalent to the GP equation 



d 1 r 



(6) 



For the case of attractive interactions, it is reasonable 
to assume that the amplitude takes a Gaussian form, 



with its width reduced due to the attractive interaction. 
The validity of this Gaussian approximation will be ex- 
amined by comparison with the numerically exact results 



in Sees. Ill and IV. Because we want to study the dy- 



namics of BEC, we allow the amplitude A(r, t) to vary 
in time as ilJTll 



A(r,t) = 



1 



TT 3 / 2 d X (t)dy(t)d Z (t) 



x exp 



2dl(t) 2dl(t) 2dl{t) 



(7) 



where di{t) (i — x,y,z) is a time-dependent real vari- 
ational parameter characterizing the width of the wave 
function along the i-axis. It can be shown that the equa- 
tion of continuity (g) and the requirement that there 
should be no mass current both at the origin and at in- 
finity uniquely determines the phase as 



(r,i) 



d x (t) 



dy(t)_^ 2 



d z (t) 



2d x (t) 2d y {ty 2d z (t) 
We thus obtain the variational wave function as |2£ 
VVar(r,i) = 



(8) 



1 



7T^d x (t)dy(t)d Z (t) 



exp 



[1 - id x {t)d x {t)} 



2 2 

^ ) [l- l dy{t)d y {t)]-^ ) 



2d%{t) 

1 - id z {t)d z {t)] 



(9) 



The imaginary terms in the exponent in Eq. (^|) describe 
the mass current associated with the change in the width 
of the wave function. To find the most probable Feynman 
path within a functional space of the complex Gaussian 
wave function (g), we substitute Eq. ([)[) into Eq. (||), 
obtaining 



Nh 



dt 



E 



x— x,y ,z 



i 



d!(t) 



d 2 {t)+di{t) 



d x {t)dy{t)d Z {t) 



(10) 



where 7 = -^=^. The stationary condition 5S/Sdi = 
gives an equation of motion for di (t) as [29| l 



di(t) = -di(t)+dr 6 (t) + 



where 



v* = \ E K 2 



7 = dVcS 

2d X (t)dy(t)d Z (t)d l {t) ^ ddt ' 

(11) 

(12) 



2d x d y d z 



may be viewed as an effective potential for the widths of 
the wave function dk(t) |l4p|. 
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III. POWER-LAW BEHAVIOR OF COLLECTIVE 
MODES 

It has been predicted in Ref. that as the criticality 
is reached, the monopole mode becomes softer according 
to the one-fourth power 1 — N/N c . In Ref. p5| , however, 
a sum rule is invoked in order to determine the numer- 
ical coefficient. Here we develop a general theory that 
determines both the power and the numerical coefficient. 
For an isotropic trapping potential, the stationary 



point of Eq. (|ll), d x 
by a positive root of 



cL, 



is determined 



(13) 



For the case of attractive interaction 7 < 0, this equation 
has a positive root if the strength of the interaction | — y | 
is smaller than the critical value 7 C MKi 



H <7c= ^74, 



(14) 



where I7I = 7 C corresponds to d st = r c = 5~ 4 / 4 . The 
stationary point d s t is always greater than r c for I7I < j c . 
The stability condition is given by 



det 



( d 2 V cS 
\ ddiddj 



= fe 4 -5)(2C t 4 + 2) 2 <0, (15) 



d=d s t 



where the derivative is evaluated at the stationary point 
d x = d y — d z = dst- It follows from Eqs. (|l5|) and Jl3| ) 
that the necessary and sufficient condition for BEC to 
have a metastable state is | — / 1 < j c . 

Let Ui(t) = di(t) — dst be the deviations of di(t) from 
its stationary point at time t. Small oscillations of these 
quantities are governed by 



Ui{t) 



E 

j=x,y,z 



d 2 V cfi 



ddiddj 



Uj(t). 



(16) 



d=d st 



Substituting Ui{t) = Ui(0)e tuJt into this yields an eigen- 
value matrix equation: 



-Ct 4 


- 3 


c t 4 - 


Ct 4 - 


1 


^-Ct 4 


Ct 4 - 


1 


Ct 4 - 



1 



to 







u x 








Uy 


= 


3 









(17) 

For this equation to have nontrivial solutions, we should 
have 



det I uj 2 Si 



•err 



ddiddj 



(co 2 



d=d at 
,-4 



5)( W 2 -2d s - 4 " 2)^ = 0. 



(18) 



Hence, we have the frequency of the monopole mode com 
and that of the doubly degenerate quadrupole mode uq 
as MM 
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FIG. 1. The effective potential f(r) = 3r 2 + 3r~ 2 + -yr~ 3 
for the width r of BEC with attractive interactions for 
| 7 | = 7c = 8 ■ 5" 5/4 , | 7 | = 0.9 7c , and | 7 | = l.l 7c . The 
dashed line shows the inflection point r c at | 7 | = 7c . 



W M = V 5 - rf s t 4 > 



^Q-V 2 ( 1 + Ct 4 )- 



(19) 
(20) 



The eigenvector corresponding to the monopole mode is 
given by (u x ,u y ,u z ) — (1,1,1), and therefore d x (t) — 
d y (t) = d z (t) = r(t) always hold. Substituting this into 
Eq. ( [Toj ) yields an effective action for the monopole mode, 

Sm = ^J dt [3f 2 - /(r) + 3f] , (21) 



where 



f{r) =3r- 2 +3r 2 +-fr- 3 (22) 



may be viewed as an effective potential for r(t). 

The forms of the function f(r) in the vicinity of | 7 | = 
7 C are illustrated in Fig. [I] [^9 22|. The potential f (r) is 
a monotonically increasing function when | 7 | > j c and 
has a local minimum when | 7 | < 7 C , indicating that j c is 
the critical strength of the interaction above which the 
condensate collapses upon itself. In other words, when 
the number of atoms exceeds the critical value 



2\/2nd s 
5 5 / 4 |a| 



(23) 



the parameter r(t) slides down the potential in time to 
r(t) — * 0, and the central density of the wave func- 
tion grows unlimitedly. In contrast, when the number 
of atoms is smaller than A^ c , f(r) has a local minimum 
at which a metastable BEC is formed. 

For I7I slightly below 7c , we set | 7 | = j c — ^7 and 
dg t = r c + 6d, and expand Eqs. ©, ©, and @ in 
small quantities ^7 and Sd. We then obtain 



3 



Sd 2 ( N\ 

(25) 
(26) 

Substituting Eq. (|24l) into Eqs. (|2j) and @, we obtain 



fid 

cj m = W20 — , 

V T r 



LU Q = 4/12-40 — . 



1/4 



\ 



J2 -s, /ID | 1--^- 



(27) 



(28) 



The result ( p7| ) agrees with that of Ref. ]l5| up to the 
numerical coefficient. 

Figure |^ compares the analytic results (|l^) and ( |20| ) 
(solid curves) with those (dashed curves) obtained by 
numerically diagonalizing the Bogoliubov equations |H| . 
Since the Gaussian approximation overestimates the crit- 
ical number of atoms N c , the numerically exact value of 
N c is used for the dashed curves in Fig. In Fig. || 
(a), \fh and y/2 indicate the Thomas-Fermi limits for the 
monopole and quadrupole modes, respectively [^o| . The 
numerical and analytic results are in excellent agreement 
for an extensive range of 1 — N/N c , and the power law 
of wm for 1 — N/N c 1 described in Eq. (^7j ) is clearly 
illustrated in Fig. || (b). 
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IV. COLLAPSING DYNAMICS OF BEC 

We assume that BEC is initially prepared in a 
metastable state with 7 slightly smaller than j c . We 
then suddenly tighten the trap potential or increase \a\ 
by a technique such as the Feshbach resonance so 
that I7I = 7 C + £7 slightly exceeds 7 C . BEC then be- 
gins to collapse. Although the collapsing behavior and 
its time scale are numerically studied in Refs. [g^ , g3| , g6| , 
no analytic forms have been reported for them. 

The effective Lagrangian for the monopole mode is ob- 
tained from Eq. ( |2l| ) as 

£M = ^(3r 2 -3r- 2 -3r 2 + | 7 |r- 3 ), (29) 



and the equation of motion reads 

r{t)=r-3{t)-r{t)-\X-\t) 



(30) 



Integrating Eq. (30) with initial conditions r(0) = r c and 
f(0) =0 gives 



dr 
di 



l[f(r c )-f{r)}. 



FIG. 2. The frequencies of the monopole and quadrupole 
modes wm and ujq, normalized by the trap frequency ujo as a 
function of N/N c sgn a, where sgn a = 1,0 and —1 for a > 0, 
a = 0, and a < 0, respectively. The solid curves show the re- 
sults of the variational method, and the dashed curves show 
the results obtained by numerically diagonalizing the Bogoli- 
ubov equations. The dotted lines at v5 and y/2 show the 
Thomas-Fermi limits of the monopole and quadrupole modes, 
respectively. The curves in (b) show the same data as in (a) 
on a logarithmic scale to clarify the power-law behavior of the 
monopole mode. 



Taking r(t) — r c — x(t) and expanding the right-hand 
side of Eq. (|3l]) up to the third power of x(t), we obtain 

( = 5<W«) + 2 ■ 5 5 /^ 7 z 2 (t) + ^l x ^(t). 
\ at J 3 

(32) 

The initial stage of the collapse, where x(t) <C 1, is dom- 
inated by the first two terms on the right-hand side of 
Eq. ( |32l ) . Integrating Eq. (|32|) by keeping only these two 
terms, we obtain 



(31) 



4 



x(t) 



4-5 



f cosh V2 • 5 5 / 4 (5 7 t - 1 
5 1 / 4 V 



-It 



(33) 



In contrast, the final stage of the collapse is dominated 



by the last term in Eq. (32) if 67 -C 1, so that we obtain 

3 ■ 5- 5 / 4 



x{t) 



(^collapse 



(34) 



where the constant of integration i C oiiapse may be inter- 
preted as the time scale for BEC to collapse. Because 
^collapse is dominated by the slow initial stage of the col- 
lapse, its evaluation requires the inclusion of the first and 
third terms of Eq. (|32|), obtaining 



t 



collapse 



d.v 



5(572; 



(3/10) 

4v^F \A 



N 

1 

Nr. 



-1/4 



1.37 



N 



1 



-1/4 



(35) 



It is interesting to note that this collapse time obeys the 
same power with respect to |1 — N/N c \ as that of for 
N just below N c (see Eq. @). 

To check the validity of these analytic results, we 
numerically integrate the time-dependent GP equa- 
tion using the finite-difference method with the Crank- 
Nicholson scheme B . Figure [| shows the time evolution 
of the peak height of the wave function \ip(r = 0,i)|, 
where we prepare an initial metastable state for I'y ] / Tc > 
1 — 10~ 7 , suddenly increase I7I to I7I/7C = 1 + 10~ 2 at 
t = 0, and then let the state evolve in time according to 
the GP equation (^) without atomic loss or to the GP 
equation ( p8|) with the atomic loss. The dashed curve is 
obtained by solving Eq. (||), and the solid curve is ob- 
tained by solving Eq. (|3^) which includes the atomic- loss 
processes due to two-body dipolar loss and three-body 
recombination loss. In the initial stage of the collapse, 
the peak density grows slowly, and at t ~ 2.52 the rapid 
implosion breaks out, where its blowup is shown in the 
inset. When the atomic loss is not included, the peak 
density increases to infinity (dashed curve). 

Figure |J shows the profiles of the wave functions 
\tp(r, t)\ obtained by numerically solving the GP equa- 
tion without loss processes (@) (solid curves) and the cor- 
responding Gaussian functions that has the same width 
(f) (dashed curves). Even at t = 0, the wave function 
deviates from the Gaussian function, particularly near 
the center of BEC. This deviation is the origin of the 17 
% discrepancy in N c between the numerically exact value 
and that obtained using the Gaussian trial wave function. 
During the initial stage of a gradual increase in the peak 
density, the deviation between the numerically obtained 
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FIG. 3. Time evolution of the peak value \ip{r = 0,t)| 
of the wave function. The solid curve shows a case with 
two-body dipolar loss and with three-body recombination 
loss (Eq. (pq)), and the dashed curve without these losses 
(Eq. (|^)). We first prepare BEC in a metastable state slightly 
below the critical point, and at t = 0, increase N so that 
N/N c — 1 = 10~ 2 . The inset shows a blow-up view of the 
intermittent implosion. 



r=2.0 




FIG. 4. The solid curves show the profiles of the wave 
functions |^)(r, £)| at t = 0,2,2.5, and 2.519, obtained by nu- 
merically solving the Gross-Pitaevskii equation (^j) for the 
same conditions as in Fig. ^. The dashed curves show Gaus- 
sian functions having the same (r) as each wave function (for 
the definition of (?), see the text). 



and Gaussian wave functions does not grow so much (see 
Fig. |] (b)). At and after the outbreak of the implosion 
(Figs. [| (c) and (d)), however, the deviation becomes 
significant, and the Gaussian approximation apparently 
breaks down, as the implosion occurs in the small region 
at the center. The implosion occurs suddenly, so we may 
define iimpiosion as the time at which the rapid implosion 
occurs. 

Until the rapid implosion takes place, the Gaussian 
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FIG. 5. (a) The time development of (f (0)) - <r(t)) ob- 
tained by numerically solving the Gross-Pitaevskii equation 
(|) for N/N c - 1 = 10 -2 , 1CT 3 , 1(T 4 , and 1CT 5 . BEC is 
initially prepared in a metastable state for N just below N c , 
and N is suddenly increased above N c . (b) The time develop- 
ment of [(f(0)} — (r(t))] _1 ^ 2 obtained by numerically solving 
the Gross-Pitaevskii equation (§) for N/N c - 1 = 1CT 3 , 1CT 4 , 
and 10 -5 , as functions of ti mp iosion — t, where ti mp iosion is the 
time at which the rapid implosion occurs. The dashed line 
shows the slope of (3 ■ 5~ 5 ^ 4 2/y/n)~ 1 ^ 2 ~ 1.49 for reference. 



trial wave function fairly well describes the collaps- 
ing dynamics. Figure ^| shows the time evolution of 
(f(0)) — (f(t)) at the initial and final stages of the collapse, 
which is obtained by numerically solving the lossless GP 
equation (||). The expectation value (f (t)) is defined as 



r# V ar(r,i)| 2 dr 



=r(t), 



(36) 



where ip vav (r,t) denotes the Gaussian variational wave 
function (^|) and r(t) is the variational parameter, and 
hence (r(0)) — (f(t)) corresponds to 2x(t)/y / iT, where 
x(t) = r c — r(t). The initial changes in the width 
of the wave functions, shown in Fig. || (a), fit the 
square of t well, and the prefactors are proportional to 
N/N c —l, in agreement with Eq. (|33|). Figure [| (b) shows 
[(f(0))-(f(t))]- 1 / 2 at the final stage of the collapse. The 
curves become linear as they approach t = iimpiosion, with 
the slope being close to (3 • 5~ 5 / 4 • 2/0F)- 1 / 2 ~ 1.49, in 
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FIG. 6. The time it takes BEC to collapse as a function 
of N/Nc — 1 for N above iV c . The open circles represent 
times at which the rapid implosions occur, all of which are 
obtained by numerically solving the Gross-Pitaevskii equation 
([]). The dashed line shows Eq. (^), and the dotted curve 
shows Eq. @ with At = 1.83. 



agreement with the inverse-square behavior in Eq. (|34|). 
The curves end at a finite value [(r(0)) — (f(t))}~ 1 / 2 ~ 2 
at t = timpiosion, since the rapid implosion occurs at this 
time. 

The time timpiosion is determined by numerically solv- 
ing the GP equation (|), while i C oiia P sc, given in Eq. (|3|), 
is analytically derived by using the Gaussian approxi- 
mation and the small-deviation approximation (|3^). Al- 
though these approximations break down just before the 
implosion, t C oiiapsc should still give the time scale of the 
collapsing dynamics, for the time it take BEC to collapse 
is dominated by the slow initial stage, where the approxi- 
mations are valid. Quantitatively, however, i co n a p Se over- 
estimates the correct implosion time ^implosion because 
Eq. (|3^) neglects higher-order terms that accelerate the 
implosion at the final stage of the collapse. Figure ^ 
shows tcoiiapso (dashed line) and numerically obtained 
implosion (open circles). Since the difference between 
implosion and £ C oiiapsc comes from the final stage of the 
collapse and we consider the case in which the initial 
numbers of atoms are almost the same, their relation 
may be given by 



^implosion — ^collapse At, 



(37) 



where At is a constant that should be determined from 
a numerical analysis. The dotted curve in Fig. |^ rep- 
resents Eq. (|37]) with At = 1.83, which is in excellent 
agreement with the numerically obtained exact iimpiosion- 
For small N/N c — 1, timpiosion asymptotically approaches 
^collapse (dashed line) , which is governed by the power law 
oc (N/N c - l)- 1 / 4 . 
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V. INTERMITTENT IMPLOSION AND 
EXPLOSION 

If the atomic loss is not included in the GP equation, 
the peak density grows to infinity, as shown in Fig. ^ 
(dashed curve). In actual experiments, when the density 
becomes very high, the atomic loss due to inelastic colli- 
sions becomes important. We, therefore, employ the GP 
equation with loss processes ]2l| as 

= + y * + <#l> - 5 ( + f 1 ^' 4 

(38) 

where L 2 and L 3 denote the two-body dipolar and three- 
body recombination loss-rate coefficients, respectively. 
The two-body (three-body) loss-rate coefficient must be 
divided by two (six) because of Bose statistics | f32f . In 
Eq. j38|), L2 and L3 are made dimensionless by multiply- 
ing N/(u)ocIq) and N 2 /(uod,Q), respectively. The atoms 
and molecules produced by inelastic collisions are as- 
sumed to escape from the trap without affecting the con- 
densate. 

The solid curve in Fig. [| shows the time evolution of 
the peak height of the wave function \i/)(r — 0, t)| fol- 
lowing Eq. fl38|), where we assume a 7 Li condensate with 
N = 1260, lu = 2tt x 144.5 Hz L 2 = 1.05 x 10~ 14 
cm 3 /s @, and L 3 = 2.6 x 10-^cm 6 /s Q. The ini- 
tial condition is the same as the lossless case (dashed 
curve). During the gradual increase in the peak density, 
a few atoms are lost, which slightly delays the time of 
implosion. When the implosion begins, the role of the 
atomic loss becomes significant. The implosion stops at 
a certain density, showing pulse-like behavior (see the 
inset of Fig. ^), and occurs several times intermittently 
in rapid sequence J27J. This behavior may be qualita- 
tively explained as follows. When the rapid implosion 
occurs and the peak density satisfies g\ip\ 2 ~ L 3 \ifj\ 4 /12, 
i.e., \ip\ 2 ~ 12g/L 3 , the collisional loss rate begins to 
surmount the accumulation rate of atoms at the cen- 
ter. When the atoms at the peak are suddenly re- 
moved due to atomic loss, the attractive force within 
a small spatial region is weakened, and the atoms be- 
gin to explode rather than implode due to the zero- 
point kinetic pressure. When the inward flow outside 
the region of implosion is sufficient to replenish the peak 
density, the subsequent implosion occurs. The intermit- 
tent implosions thus occur as the result of a competi- 
tion between the loss of atoms and their accumulation 
due to the attractive interaction; therefore, these implo- 
sions should be distinguished from other oscillatory be- 
haviors due to various mechanisms such as collapse-and- 
growth cycles |2^,^,^l| and the small collapses described 
in Ref. pl| . In the case of 7 Li, the peak density is esti- 
mated to be \ip\ 2 ~ 12g/L 3 ~ 400 2 , which qualitatively 
agrees with the solid curve in Fig. [| We should note 



here that the peak density |V>| 2 ~ 48Trha / (mL 3 ) , where 
the dimension is restored, is independent of the number 
of atoms N and the trap frequency ljq, indicating that 
the implosions are a local phenomenon depending only 
on the nature of the atoms themselves. 

The explosion following an implosion of atoms is pre- 
dicted in Refs. ]2l],^7j and has been experimentally ob- 
served as a broadly scattered "hot" atom cloud of about 
100 11K S. From the above discussion, the energy scale 
of the explosion might be estimated as g\ip\ 2 ~ 12g 2 /L 3 . 
jj This energy is, however, the highest one that only a small 
number of atoms at the center of BEC may acquire. In 
fact, the atoms scattered by the explosion have a broad 
momentum distribution. Figure Q (a) shows a normalized 
energy distribution of atoms scattered by the explosion, 
where the parameters of 7 Li are used; the lower panel of 
Fig. (a) shows the same distribution with L 3 is mul- 
tiplied by 100. These energy distributions are obtained 
from the wave functions at t = ^implosion + tt/2, at which 
the ejected part of BEC spreads maximally and are there- 
fore easy to be discerned from the remnant part of BEC 
that remains at the center. The energy per atom at po- 
sition r is given by 



— ) V(r)/|^(r)| 2 , (39) 



where the interaction energy is omitted since the density 
of the widespread atoms are low enough. In obtaining 
the energy distributions and the mean energies, the rem- 
nant part of the condensate at the center is excluded. 
While the energy scale is different, the two distributions 
in Fig. [t| (a) have similar shapes. The number of atoms 
scattered by the explosion is 20 ~ 40 % of the total num- 
ber of atoms. The dependence of this fraction on the 
parameters g and L3 is complicated, since it depends on 
the number of times and the detailed structure of the 
intermittent implosions that are governed by highly non- 
linear dynamics. 

Figure (b) shows the mean energies of the atoms 
scattered by the explosion as a function of g 2 /L 3 . The 
open circles are obtained using the value of g of Li with 
u>o = 2ir x 144.5 Hz and the three-body recombination 
rates L 3 , IOL3, and IOOL3. In order to see the depen- 
dence on g, the value of g is halved in the solid circles. 
The open and solid circles have almost the same depen- 
dence on g 2 /L 3 , suggesting that the explosion energy is 
determined by the ratio g 2 / L3, and not by g and L 3 sepa- 
rately. The mean energy is roughly given by ~ Q.\g 2 / L 3 , 
or restoring the dimension by 



0.1 x m-K 2 



(40) 



We note that this energy is independent of the number 
of atoms N and the trap frequency ujq, as with the peak 
density discussed above. These results indicate that the 
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FIG. 7. (a) Normalized distributions of the energy of an 
atom scattered by the explosion. The same parameters and 
the initial conditions as those in the solid curve in Fig. ^| are 
used in the upper distribution, and the three-body recombina- 
tion rate is multiplied by 100 in the lower one. (b) The mean 
energies of an scattered atom. The rightmost open circle is 
obtained using the parameters of 7 Li, and the other open cir- 
cles are obtained by increasing the three-body recombination 
rate by 10 and 100. In the solid circles, the s-wave scattering 
length is halved. 



tal result of the explosion energy ~ 100 nK based on 
Eq. (|io|). Assuming a = — 1 nm, we obtain ~ 6x 10~ 28 
cm 6 /s, which is consistent with the value of L 3 around 
170 G |3!|. The above-described method thus provides an 
independent way to determine the three-body recombi- 
nation rate in addition to the current methods ||3^.|37|,|35| . 



VI. CONCLUSIONS 

We have analytically and numerically studied the dy- 
namics of BEC with attractive interactions below and 
above the critical number of atoms 7Y C . We have shown 
that our analytic approach successfully describes the col- 
lective oscillations below N c as well as the collapsing dy- 
namics of BEC above N c until the rapid implosion takes 
place. We have found that the time scales of the dy- 
namics follow the same power law on either side of the 
critical point. Just below N c , the collective frequency 
of the monopole mode is proportional to (1 — N/Nc) 1 / 4 , 
and just above N c , the time it takes BEC to collapse is 
proportional to (N/N c — l)" 1 / 4 . 

When the atomic loss due to inelastic collisions is taken 
into account, the implosion occurs not only once but sev- 
eral times intermittently, and the implosion is then con- 
verted to an explosion. We have found that the mean 
energy of an atom scattered by the explosion is given by 
Eq. (p£)|), which is independent of the number of atoms 
and the trap frequency 
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global nature of the system is unimportant with regard to 
implosion and explosion, and that only the local nature 
of the atoms determines their behavior. Using Eq. (40), 



the mean kinetic energy of the atoms scattered by the 
explosion in the experiment at Rice is estimated to 
be ~ 80/iK. In the experiment at JILA ||, on the other 
hand, the three-body recombination rate of 85 Rb is un- 
certain in the vicinity of the Feshbach resonance, where it 
significantly depends on the applied magnetic field p5[ . 
The experiment of collapse is carried out by applying 
magnetic field of about 170 G at which sign of the inter- 
action changes. Because of uncertainty in the three-body 
recombination rate, we estimate it from the experimen- 
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